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Abstract 

Estimates  of  the  power  spectral  density  have  been  found  to  be  very  useful  in  a  variety  of  signal 
processing  applications  over  the  last  several  decades  [1,  2,  3].  Higher  order  spectra  contain  information 
not  present  in  the  power  spectrum  and  recently,  estimates  of  higher  order  spectra  have  been  shown 
to  be  useful  in  certain  signal  processing  problems  [4,  5,  6,  7,  8,  9,  10].  In  particular,  estimates  of  the 
bispectrum  and  bicoherence  have  been  found  useful  in  detecting  non-Gaussianity  and  non-linearity  [4, 
6,  9],  in  system  identification  [11],  and  in  detecting  transient  signals  [10,  12]. 

This  paper  is  concerned  with  the  estimation  of  the  bispectrum  and  bicoherence  of  underwater 
acoustic  signals.  We  first  introduce  and  describe  some  of  the  properties  of  the  bispectrum  and  bico¬ 
herence.  The  bispectral  estimation  algorithm  as  implemented  here  at  the  Marine  Physical  Laboratory 
is  then  described  in  detail.  Finally,  several  bispectra  are  estimated  from  actual  data  and  the  results 
are  analyzed. 

In  the  results  portion  of  the  paper  we  concentrate  on  three  particular  properties  of  bispectrum 
estimation.  We  show  how  the  bispectrum  estimate  can  be  used  to  detect  non-Gaussianity,  non-linearity, 
and  harmonic  coupling.  The  detection  of  harmonic  coupling  is  shown  to  be  a  useful  property  of  the 
particular  bispectrum  estimation  algorithm  employed  here. 
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1  Introduction 


Power  spectrum  estimation  has  been  a  useful  signal  processing  tool  for  many  years  [1,  2,  3].  Since  the 
power  spectrum  is  related  to  the  autocorrelation  it  does  not  contain  any  information  about  moments  higher 
than  the  second  order.  It  is  therefore  of  limited  utility  when  one  is  trying  to  characterize  or  detect  non- 
Gaussian  noise.  In  system  identification  problems  power  spectral  estimation  techniques  can  only  determine 
the  magnitude  response  of  the  system.  The  phase  response  is  unknown.  (It  is  often  assumed  however  that 
system  is  minimum-phase.  If  this  assumption  is  correct  then  the  complete  response  can  be  determined 
using  power  spectrum  estimation  methods.) 

The  bispectrum1  is  a  third-order  spectrum  which  has  a  number  of  properties  that  potentially  make 
it  a  valuable  signal  processing  tool  [13,  14,  15].  In  practice  the  bispectrum  has  proven  to  be  particularly 
useful  in  detecting  non-Gaussianity  and  nonlinearity  of  ambient  noise  [4,  6,  7].  More  recent  work  [10,  12] 
indicates  that  the  bispectrum  can  also  be  useful  in  transient  detection.  The  bispectrum  can  also  be  used 
to  determine  the  phase  response  in  system  identification  problems  [11], 

The  remainder  of  this  section  is  used  to  introduce  notation  and  to  define  the  bispectrum  and  bicoher¬ 
ence.  In  Section  2  we  describe  in  detail  the  particular  bispectral  estimation  algorithm  implemented  here 
at  the  Marine  Physical  Laboratory.  Certain  properties  of  the  bispectrum  and  bicoherence  are  described  in 
Section  3.  We  discuss  a  particular  property  of  our  bispectrum  estimator  which  is  useful  in  determining  if 
spectral  lines  are  harmonically  related.  We  also  discuss  some  of  the  statistical  properties  of  the  bispectrum. 
These  properties  are  extremely  useful  in  interpreting  bispectral  estimation  results.  In  Section  4  we  present 
results  of  bispectral  estimation  analysis  as  applied  to  data  collected  from  a  Swallow  float  and  also  to  data 
collected  from  single  elements  of  the  MDA  array. 

1.1  Definitions 

Let  x(t)  represent  a  real  zero-mean  stationary  random  process. 

1.1.1  Moments 

The  first-order  moment  (or  mean)  is  defined  as 

mx  -  £'{x(i)},  (1) 

which,  for  the  stated  assumptions,  is  constant  and  equal  to  zero.  The  second-order  moment  (or  autocor¬ 
relation)  is 

r**(r)  =  E{x(t)x(t  +  r)}.  (2) 

The  bispectrum  is  related  to  the  third-order  moment  (or  autobicorrelation)  which  is  defined  as 

Rxxx(n,  r2)  =  E{x(t)x(t  +  Ti)x(t  +  t2)}.  (3) 

For  the  discrete  sequence  x[n]  obtained  by  sampling  the  continuous-time  signal  x(t)  every  T  seconds,  i.e. 


a:[n]  =  x(nT),  the  corresponding  definitions  are 

mx  = 

£{*[«]} 

(4) 

r*®[n]  = 

£{a;[m]zr[m  +  n]} 

(5) 

RXxx[m,  n]  = 

£{a:[/]x[/  +  m]x[l  +  n]}. 

(6) 

'The  function  commonly  called  the  bispectrum  is  more  properly  termed  the  bispectral  density. 
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1.1.2  The  Power  Density  Spectrum  and  the  Bispectrum 

The  power  density  spectrum  is  defined  as  the  Fourier  transform  of  the  autocorrelation  function  i.e., 

/oo 

rxx(r)e~j2wfT  dr.  (7) 

■CO 

Assuming  that  x(t)  is  a  voltage  then  the  units  of  Pxx{f)  are  Volts2/Hertz  (or  assuming  the  power  is 
dissipated  in  a  1  Ohm  resistance,  Watts/Hertz). 

The  bispectrum  is  defined  as  the  Fourier  transform  of  the  autobicorrelation  sequence, 

/OO  [OO 

/  Rxxx(r1,T2)e-^^+^drldr2  (8) 

■CO  J  —  OO 

with  corresponding  units  of  Volts3/Hertz2. 

For  discrete  sequences  the  corresponding  definitions  are 

PM  =  T  £;  rxx[n]e-^nT  (9) 

n=-oo 

CO  00 

JW/1./2)  =  E  E  fl***[™>«K;'2’(/,m+/3n)T-  (10) 


The  value  of  the  bispectrum  depends  on  the  amplitude  of  the  component  signals,  and  it  is  desirable  to 
have  a  normalized  form  of  the  bispectrum  [8].  The  bicoherence  is  defined  as  [8] 


bXXX  (fuh) 


_ Bxxxjfl  1  f 2) _ 

\/PM)PM)Pz*(fl+f2)‘ 


(11) 


The  bicoherence  has  units  of  Hertz-1/2 


* 
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2  Estimation  of  the  Bispectrum 

Bispectral  (and  power  spectral)  estimation  techniques  can  be  classified  as  being  either  parametric  or  non- 
parametric  (conventional)  and  also  as  being  either  direct  or  indirect  [14].  Parametric  estimators  use  a 
model  (usually  of  the  AR,  MA,  or  ARMA  type)  for  the  process  under  investigation  and  obtain  an  estimate 
of  the  bispectrum  by  first  estimating  the  parameters  of  the  model.  Non-parametric  estimators  make  no 
assumptions  about  a  model  for  the  process.  Parametric  techniques  can  provide  better  estimates  of  the 
bispectrum  than  non-parametric  techniques  in  situations  where  an  assumption  about  the  underlying  model 
is  valid  [15]  In  particular,  parametric  estimation  techniques  are  often  used  in  the  system  identification 
problem.  Direct  type  algorithms  estimate  the  bispectrum  directly  from  the  data,  indirect  types  compute 
the  bispectrum  from  an  estimate  of  the  autobicorrelation. 

The  particular  algorithm  presented  here  falls  into  the  non-parametric,  direct  estimation  class  of  bis¬ 
pectrum  estimators.  It  has  many  similarities  to  Welch’s  method  of  power  spectrum  estimation,  and  so  a 
brief  overview  of  that  technique  is  presented  after  a  discussion  of  notation. 

2.1  Notation 

Let  x[n\  denote  the  nth  sample  of  the  random  process  x(t),  i.e. 

a?[n]  =  x(nT)  0  <  n  <  N  —  1,  (12) 

where  T  is  the  sampling  interval  and  it  is  assumed  that  we  have  N  data  points.  In  order  to  reduce  the 
variance  of  a  spectral  estimate  the  data  is  usually  divided  into  segments  (unfortunately  this  also  reduces 
the  resolution).  Assume  that  the  data  is  divided  into  P  segments  of  D  samples  with  a  shift  of  S  samples 
between  segments.  The  weighted  pth  segment  can  be  written  as 

x(p)[n]  =  tu[n]x[n  +  p5]  0<n<D-l.  (13) 

The  discrete-time  Fourier  transform  (DTFT)  of  the  weighted  pth  segment  of  data  is  equal  to 


o-i 

x«(,)=r£  x^[n)e~j2vfnT.  (14) 

n=0 

The  factor  of  T  is  included  in  the  definition  so  that  the  discrete-time  Fourier  transform  will  approximately 
equal  the  corresponding  continuous-time  transform.  The  discrete  Fourier  transform  (DFT)  of  the  same 
segment  is  equal  to 

o-i 

X<*>[*]  =  Y,  x^[n}e-^Tnk'D.  (15) 

n=G 

The  factor  of  T  is  omitted  in  the  definition  of  the  DFT  to  conform  to  standard  notation. 


2.2  Welch’s  Method  of  Power  Spectrum  Estimation 

There  are  several  methods  for  estimating  the  power  density  spectrum.  The  various  methods  are  discussed 
in  detail  in  several  texts  on  signal  processing  [1,  2,  3].  The  bispectrum  estimation  procedure  presented  in 
the  next  section  is  analogous  to  Welch’s  method  for  power  density  spectrum  estimation  and  so  the  equations 
which  define  that  method  will  be  presented  here.  The  notation  closely  follows  that  of  Marple  [1]. 

The  sample  spectrum  of  the  weighted  pth  segment  is  given  by 

HPJ(f)  =  ^l*(p)(/)|2,  (16) 
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where 


(17) 


D-l 


Up  —  T  ^2,  U,2W- 


n=0 


The  Welch  estimate  is  obtained  by  averaging  the  sample  spectrum  from  each  segment 


p-i 


PrM)  =4E 


P—0 


If  a  DFT  is  used  to  estimate  the  spectrum,  the  equations  become 

HPM  =  §V(P)MI2 


p-i 


p*m  = 


(18) 

(19) 

(20) 


P—0 


where  X^p^[ifc]  is  the  DFT  of  the  pth  weighted  data  segment. 

For  sinusoidal  signals  it  is  often  more  convenient  to  work  with  a  function  that  indicates  the  power  in 
each  DFT  bin.  Denoting  this  function  by  Vxx[k]  we  have 


vim  = 


p 

p~  i 


vm  = 


(21) 

(22) 


p= o 

Vxx[k]  has  units  of  Volts2  (or  Watts  if  measured  across  a  1  Ohm  resistor). 

2.3  Bispectrum  Estimation 

Our  bispectrum  estimation  algorithm  is  based  upon  the  following  equation 

1 


Bxxx{fi,fi)  =  lim  E{ 


where 


M—*oo  1  (2  M  +  1  )T 
M 


XM(f)  =  T  £  x[m]e-^nT. 


(23) 

(24) 


m——M 


The  proof  of  this  equation  is  fairly  straightforward.  After  using  the  definition  of  X\i(f)  in  Eq.23  we  have 

M  M  M 


Bxxx{fuh)  =  Jim 


T2 


M—co  (2 M  +  1)  _ 


Y  Y  Y  Rxxx[m-l,n-l\e-^^m-l^Te-^^n-^T.  (25) 


m=— M  n=— JV/ /=— M 


It  can  be  shown  that 

M  M  M  2  M  2  M 

X  Y  XI  /[m -/,«-/]=  X  X  (2Af+  1-S[m,n])f[m,n],  (26) 


mzz-M  n  =  -M  l=-M 


m=— 2M  n  —  —  2M 


where 


S[m,  n]  =  < 


2  M  +  1 

m  >  0;  —2 M  <  n  <  —2 M  +  m 

m  —  n 

m  >  0;  —2 M  +  m  <  n  <  0 

m 

m  >  0;  0  <n  <m 

n 

m  >  0;  m  <  n  <  2 M. 

(27) 
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The  values  of  S[m,  n]  for  m  <  0  can  be  obtained  from  S[— m,  — n]  =  S[m,  n].  Note  that  S[m ,  n]  is  positive 
and  is  less  than  or  equal  to  2 M  +  1.  Using  this  in  Eq.  25  gives 

Bxxxifl  1/2)  = 


where  it  has  been  assumed  that 

CO  00 

T.  T  S[m,  n]Rxxx[m,  n ]  <  00.  (31) 

m=— 00  n  =  — 00 

The  last  equation  for  the  bispectrum  is  identical  to  Eq.  10  and  therefore  proves  Eq.  23. 

This  proof  suggests  that  the  following  algorithm  for  estimating  the  bispectrum.  Let 

1  P_1 

JW/i -  h)  =  p  £  #&(/i  -  /a)  (32) 

p=0 

where  BxPxx(fi ,  h)  is  the  sample  bispectrum  of  the  weighted  pth  data  segment 

B'ilifuh)  =  ~X^\h)X(p\f2)X^{h  +  f2)  (33) 

and  X(p^(f  )  is  the  discrete-time  Fourier  transform  as  defined  previously  in  Eq.  14.  In  Eq.  23  t4  is  equal 
to  1/(2 M  +  1)T.  It  is  difficult  to  say  what  t/j  should  be  now  because  we  are  using  windowed  data  that  is 
finite  in  extent  and  our  bispectrum  estimate  may  be  biased.  The  remainder  of  this  section  is  concerned 
with  finding  the  value  of  t/j  which  gives  an  unbiased  estimator. 

Assuming  that  E{BipJx{fi,  /2)}  does  not  depend  on  p  (this  is  verified  later)  then  the  mean  value  of 


BXxx(fi,h)  is  equal  to 

E{Bxxx(fuh)}  =  E{BipJx(fuf2)}.  (34) 

Substituting  Eq.  33  yields 

E{BXxx{h,h)}  =  ~E{X^{h)X^p\h)XW9{h  +  /2)}  (35) 

or  (after  using  Eq.  14) 

E{Bxxx{fi,h)}  =  7T  E  E  E  E{x(P)[/]x(P)Hx^)[n]}e^'2,r^(m-0+A(n-01T-  (36) 

^  /=0  m= 0  n=0 

The  expectation  on  the  right  side  of  this  equation  is  equal  to 

£{x^[/]x^[m]x^[n]}  =  RXXx{™  —  l,n  —  /]u>[/]u;[m]u;[n].  (37) 

Substituting  into  Eq.  36  gives 

E{Bxxx{h,  h)}  =  Tf  (38) 

6  ;=0  m=0  n=0 


T 2 


ZIYI  ZIYI 


lim  ,  . .  . 

M-00  (2 M  +  1)  £ 


y  y  (2M  +  1  -  S[m,  n])J?c„[m,  n]e~^^m+^T  (28) 


m=-2M  n  =  — 2A/ 

2  M  2  M  cr  ■« 

2  V'  /1  o[m,  n] 


lim  T2  y  y  (1- 

A/—00  ^  v 


2M  +  1 

T2  y  y  Rxxxlm^le-WBm+h^T 


)Rxxx[m,  n]e-i2*V'">+hn)T 


m=-2Af  n=— 2A/ 
00  00 


(29) 

(30) 


m=— co  n=  — 00 
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The  discrete-time  third-moment  function  is  the  inverse  Fourier  transform  of  the  discrete-time  bispectrum 

pl/(2T)  ,1/(2  T) 


Rx 


/  Bxxx(Xu  dXx  dA2. 

■1/(2 T)  J  — 1/(2T) 


(39) 


Substituting  into  Eq.  38  and  simplifying  gives  (the  limits  of  integration  are  omitted,  but  understood  to  be 
from  —1/(2 T)  to  1/(2 T)  in  the  remainder) 

E{Bxxx(f1,f2)}  =  ^J  j  Bxxx(\1,X2)W(f1-\1)W(h-\2)W'(f1+f2-\l-\2)d\1d\2  (40) 

0—1 

W(f)  =  T^2  w[n]e~j2rfnT ■  (41) 


where 


n=  0 


This  can  be  written  more  clearly  as 

E{Bxxx(fi,h)}  =  j-J  j  Bx„(\u\2)*(fi  -A1,/2-A2)dAi  d\2 


(42) 


where 


*(/i,/a)  =  W{h)W{h)W*{h  +  /*)•  (43) 

The  last  equation  for  E{Bxxx(fx,  f2)}  indicates  that  the  estimate  is  proportional  to  the  true  bispectrum 
convolved  with  $(/i,/2),  i.e., 


E{Bxxx(fuh)}  =  Jj-BxxxUuh)  *  *(/l,/2), 


(44) 


and  so  our  estimate  is  a  biased  estimator  of  the  true  bispectrum. 

It  is  not  possible  to  find  a  value  of  Ub  that  gives  an  unbiased  estimate  of  the  bispectrum  for  all  frequency 
pairs.  However,  it  is  possible  to  find  a  value  which  gives  an  unbiased  estimate  of  the  mean  cubed  process. 
The  third-moment  sequence  corresponding  to  our  bispectrum  estimate  is 


E{Rxxx[m,n]}  =  -^-Rxxx[m,n]'H[m,n] 
Ub 


(45) 


where  4([m,  n]  is  the  inverse  Fourier  transform  of  4>(/i ,  f2).  By  setting  Ub  equal  to  '£[0,0]  an  unbiased 
estimate  of  #rrX[0,0]  (the  mean  of  the  cubed  process)  is  obtained.  Ub  (^[0,0])  can  be  determined  from 


Ub  =  /  j$(fufi)dfidf2 

=  J  J  W(fx)W(f2)W*(fx  +  h)  dfx  df2 

=  T3  £  £  £  /  e-i2*h{rn-i)T  dfx  f  e-i**h(n-l)T  df2 

I  ^ - n - n  j 


1=0  m=Q  n= 0 
0-1 


=  t£u,3[(], 


(46) 

(47) 

(48) 

(49) 


1=0 

the  last  expression  being  the  most  convenient. 

If  a  DFT  (see  Eq.  15)  is  used  to  estimate  the  bispectrum,  the  corresponding  equations  are 

B{£x[hM  =  ^-X^[kx}X^{k2}X^[kx  +  k2] 

Ub 

P-1 


Bxxx[kx,k2]  = 


(50) 

(51) 


p= o 
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The  three  previous  equations  provide  the  basis  for  BISPCT,  a  computer  program  which  estimates  the 
bispectrum  of  a  process  from  its  sampled  data  [16]. 

When  working  with  sinusoidal  signals  it  is  often  convenient  to  define  a  scaled  version  of  the  estimated 
bispectrum  that  is  equal  to  (see  Eq.  22) 

Bxcx[ki,  £2]  =  ^jyp^2  ^2] •  (52) 

From  the  preceding  equations  we  have 


1  P_1 

Bxxx[kuh]  = 


(53) 

(54) 


The  units  of  Bxxx[ki,  Jfe2]  are  Volts3. 


2.4  The  Bicoherence 

We  can  estimate  the  bicoherence  using  estimates  of  the  bispectrum  and  the  power  density  spectrum, 

(55) 

The  magnitude  of  the  bicoherence  is  a  convenient  test  statistic  for  use  in  the  detection  of  coupled 
processes. 


bXxx(fi )  f*i)  — 


_ Bxxxjfl  I  ft) _ 

f  P**(h)P*x(h)Pxx(fl  +  h) 
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3  Properties  of  the  Bispectrum  and  Bicoherence 


3.1  Input/Output  Relationships 


If  x(t)  is  the  input  to  a  linear  system  with  impulse  response  h(t)  and  y(t)  is  the  corresponding  output, 
then  the  third-order  moment  of  y(t)  is  related  to  that  of  x(t)  by  [11] 


Ryyy(Tl  >  ^2) 


-IS! 


Rxxx{r 1  +  t3  -  ti,Ti  +  t3  -  t3)h(ti)h(t2)h(t3)  dti  dt2<it3. 


If  we  let  H(f)  represent  the  transfer  function,  the  corresponding  bispectra  are  related  by  [11] 

Byyy{h,h)  =  Bxxx{h  ,  f3)H  (/l)tf  (/2)iT ’(/t  +  /2)- 
The  relation  between  the  bicoherence  functions  is 


byyyifu  /z)  —  bxxx(fh  /z) 


H{h)  H(h)  £Ui±M 
\H(fl)\\H(f2)\\H(f1  +  f2)\- 


(56) 

(57) 

(58) 


Note  that  the  output  bicoherence  is  affected  only  by  the  phase  of  the  transfer  function.  (The  magnitude 
of  the  bicoherence  of  the  output  is  equal  to  the  magnitude  of  the  bicoherence  of  the  input.) 


3.2  Symmetries  in  the  Bispectrum  Plane 

The  third-order  moment  sequence  of  a  stationary  process  has  the  following  properties  [11] 

Rxxx(t1i  T2)  —  Rxxx{T2i  ri)  =  Rxxx(.~T 2,  ~  7"z)  =  Rxxx(~ fl ,  ~Tl  +  T2) 

=  Rxxx(-n  +  r2,  -n)  =  RXXx(r 1  -  r2,  — r2).  (59) 

If  the  function  Rxxx(ti  ,  r2)  is  known  in  any  of  the  six  regions  of  Fig.  1  this  set  of  equalities  can  be  used  to 
determine  it  everywhere.  From  the  above  symmetry  relations  and  the  definition  of  the  bispectrum  (Eq.  8) 
it  follows  that 


BXxx(fi,  h)  —  Bxxx(f2,fi)  —  Bxxx(  fi  —  /2>/i)  —  Bxxx(—fi  —  f2,fi) 

—  Bxxx(f 2,  fi  /2)  =  Bxxx(fi,  —  fi  —  /2). 

From  Eq.  8  it  also  follows  that  the  bispectrum  obeys  the  following  conjugate  symmetry  relation 

Bxxx(fl,f2)  =  B*xxx(-f\,-f2)- 


(60) 

(61) 


These  relations  imply  that  if  the  bispectrum  is  known  in  any  of  the  twelve  regions  shown  in  Fig.  3.2  it  can 
be  determined  everywhere  in  the  plane. 

From  Eq.  10  it  can  be  seen  that  the  bispectrum  of  a  discrete  sequence  is  doubly  periodic  with  period 
1/T: 

Bxxx(fi,  /2)  =  Bxxx(fi  +  m/T,  /2  +  n/T).  (62) 


3.3  Bispectra  of  Gaussian  Processes 

The  autobicorrelation  of  a  stationary,  Gaussian  process  is  equal  to  zero  which  implies  that  the  correspond¬ 
ing  bispectrum  is  also  zero.  This  property  allows  the  bispectrum  to  be  used  as  a  tool  for  the  detection 
of  non-Gaussian  processes.  Note  that  a  zero  bispectrum  does  not  imply  a  Gaussian  process,  i.e.  processes 
which  are  not  Gaussian  may  have  zero  bispectra.  For  example,  any  process  whose  samples  are  statistically 
independent  and  have  a  symmetric  density  function  will  have  a  zero  autobicorrelation  sequence  and  hence 
a  zero  bispectrum. 

Hinich  [17]  has  developed  a  test,  based  on  the  magnitude  of  the  bicoherence  estimate,  for  detecting 
non-Gaussianity. 
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3.4  Bicoherence  of  Linear  Processes 

A  linear  process  is  a  time  series  that  can  be  expressed  in  the  form  [6] 

(63) 


00 

x[n ]  =  ^  a(m)e(n  —  m) 

m=  —  oo 


where  e(n)  is  a  purely  random  series,  i.e.  all  of  the  {c(l),  e(2), . . .}  are  mutually  independent.  Theoretically, 
the  bicoherence  of  a  linear  time  series  is  constant  [7,  6]. 

A  non-zero  bispectrum  (or  bicoherence)  implies  non-Gaussianity.  A  non-constant  bicoherence  implies 
nonlinearity.  It  should  be  noted  that  this  implies  that  all  nonlinear  processes  are  non-Gaussian. 

Hinich  [17]  has  developed  a  test,  based  on  the  spectral  flatness  of  the  bicoherence  estimate,  for  detecting 
non-linearity. 


3.5  Bispectra  of  Harmonic  and  Quadratically  Coupled  Processes 

Assume  that  x[n]  consists  of  a  sum  of  three  sinusoidal  components, 

3 

x[n]  =  am  cos(2 nfmTn  +  (64) 

m=l 

where  /3  =  f\  -f  /2.  The  three  component  sinusoids  could  also  be  harmonically  related  (in  which  case 
all  three  frequency  components  will  be  multiples  of  a  fundamental  frequency).  One  situation  in  which 
the  three  components  satisfy  /3  =  /j  +  /2  and  are  not  necessarily  harmonically  related  occurs  in  what 
is  known  as  quadratic  coupling.  Quadratic  coupling  can  occur  when  two  time  series  (e.g.  sinusoids  with 
frequencies  f\  and  /2)  pass  through  a  transfer  function  that  is  nonlinear.  Assuming  we  can  approximate  the 
input /output  relationship  by  the  first  two  terms  in  a  Taylor  series  expansion  as  y(x)  =  ax  +  bx~  then,  if  the 
input  consists  of  two  sinusoids,  the  output  will  have  frequency  components  at  the  original  two  frequencies 
and  also  at  the  sum  and  difference  frequencies.  The  sum  component  at  frequency  /3  satisfies  the  above 
relationship.  In  addition,  the  phases  satisfy  the  relation  <t>3  —  <f>i  +  <j> 2-  This  is  termed  quadratic  phase 
coupling  because  it  occurs  whenever  the  output  and  input  are  related  through  a  quadratic  nonlinearity. 

The  discrete-time  Fourier  transform  of  the  uniformly  weighted  pth  segment  of  data,  evaluated  at 
frequency  fi  is  approximately  equal  to 

X(p)(/i)  =  a1DTej('t,l+2*hpTS)  (65) 

where  S  is  the  number  of  samples  that  we  shift  between  overlapping  segments.  The  corresponding  expres¬ 
sions  for  X^p\f3)  and  X00(/3)  are 

A<p)(/2)  =  a2DTeH,h+2*f*pTS'>  (66) 

X(p)(/3)  =  a3DTej<'<t‘3+2,rf3pTS\  (67) 


The  estimated  power  density  spectrum  evaluated  at  these  same  frequencies  is  equal  to  (see  Eqs.  16  and  18), 


PM)  =  ai(DT)  (68) 

Pxx(h)  =  4  (DT)  (69) 

Pxx(h)  =  a2(DT)  (70) 

(71) 

The  sample  bispectrum  (see  Eq.  33)  evaluated  at  ( f\  ,/2)  is  equal  to 

BxJx(fuh)  =  aia2a3(jDT)2e-'(?>1+<,!’3_^+2,r(/l+/2-/3)pT5).  (72) 
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(73) 


For  harmonic  and  quadratically  coupled  processes  in  which  fa  =  fa  +  fa,  this  reduces  to 

B(xpJx(fa,fa)  =  aia,MDT)2e^+^-^\ 

Note  that  the  sample  bispectrum  in  this  case  is  independent  of  the  segment  index  p. 

If  the  three  frequency  components  are  a  subset  of  lines  from  a  harmonic  process,  it  is  usually  assumed 
that  the  three  phases  are  independent  and  identically  distributed  random  variables  uniformly  distributed 
over  a  2jt  interval.  In  this  case  the  ensemble  average  of  the  sample  bispectrum  is  zero  and  therefore 
the  bispectrum  is  also  theoretically  zero  [13].  Fourth-order  moments  are  then  required  to  determine  if 
the  components  are  harmonically  related  [18].  However,  our  bispectrum  estimator  is  based  upon  a  time 
average  of  the  sample  bispectra  instead  of  an  ensemble  average.  Since  the  sample  bispectrum  does  not 
depend  upon  the  segment  (time)  index  p,  the  value  of  our  bispectrum  estimator  is 

Bxxx(fa,fa)  =  (74) 

The  paper  by  Huber  [8]  recognizes  the  difference  between  the  value  of  the  estimator  and  the  value  of  the 
true  bispectrum  in  the  case  of  harmonically  related  lines.  They  use  the  non-zero  value  of  the  bispectrum 
estimator  as  an  aid  in  determining  whether  lines  in  a  power  density  spectrum  are  harmonically  related. 
The  estimated  bicoherence  in  this  case  is  equal  to  (see  Eq.  11) 

bxxx(fi,fa)  =  VDTei<-*'+,h-*>\  (75) 

It  should  be  noted  that  in  this  particular  problem  our  algorithm  is  a  biased  estimator  of  the  true  bispectrum 
at  the  point  (fa,  fa).  Fortunately,  this  bias  can  be  used  to  detect  the  presence  of  harmonics. 

In  the  quadratic  phase  coupling  problem  the  phases  (j>\  and  fa  are  statistically  independent  and  fa  = 
<j> i  +  fa.  The  phase  of  the  sample  bispectrum  is  equal  to  zero  in  this  case  and  the  ensemble  and  time 
averages  of  the  sample  bispectra  are  equal.  The  value  of  the  bispectrum  estimator  in  this  case  is 

BXXX  (fa,  fa)  =  axdva-ziDT)2 ,  (76) 

and  the  bicoherence  estimate  is  equal  to 

bxxx(fa,fa)  =  y/DT.  (77) 

Points  at  which  the  magnitude  of  the  bicoherence  is  equal  to  \fUf  indicate  the  presence  of  either  a 
harmonic  or  a  quadratically  coupled  process. 

This  derivation  assumes  that  the  component  frequencies  are  stable  over  the  entire  observation  interval. 
If  this  is  not  the  case,  the  bispectral  and  bicoherence  peaks  will  be  reduced  in  amplitude.  Noise  will  also 
reduce  the  magnitude  of  the  bicoherence. 

3.6  Statistical  Properties  of  Bispectra 

As  we  have  seen,  a  non-zero  value  in  the  bispectrum  frequency  plane  is  an  indication  of  non-Gaussianity 
or  frequency  coupling.  Of  course,  the  presence  of  noise  will  also  create  peaks  in  the  bispectral  plane  and 
so  we  must  question  whether  or  not  a  peak  is  statistically  significant,  i.e.  indicative  of  non-Gaussianity  or 
frequency  coupling. 

Hinich  [17]  provides  a  fairly  complete  statistical  description  of  the  bispectrum  and  bicoherence  and 
describes  in  detail  tests  for  non-Gaussianity  and  nonlinearity. 

Huber,  et  al.  [8],  provide  a  simple  approximation  that  is  convenient  for  first-cut  data  analysis.  Since 
bispectral  estimates  are  asymptotically  unbiased  and  asymptotically  normal  under  mild  conditions  [19], 
they  indicate  that  the  variance  of  the  off-diagonal  elements  of  a  bispectrum  estimator  of  the  type  discussed 
here  is  equal  to 

D2T 

var{B*xr(/i,  fa)}  =  -jj-Pxx(fa)Pxx(fa)Pxx(fa  +  fa),  (78) 
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Figure  4:  Variance  of  the  bispectrum  of  white  Gaussian  noise  as  a  function  of  segment  length  (rectangular 
window).  The  top  line  and  corresponding  points  (diamonds)  are  the  theoretical  and  simulation  results 
when  no  overlapping  was  used.  The  middle  line  and  middle  set  of  points  (pluses)  correspond  to  50% 
overlap  and  the  bottom  line  and  points  (squares)  correspond  to  75%  overlap. 


where  Pxx(f)  is  the  power  density  spectrum  of  the  corresponding  process.  (The  variance  of  a  complex 
quantity  2  is  defined  here  as  var{z}  =  E{\z |2}  —  |i?{z}|2.)  This  result  is  applicable  only  for  the  case  in 
which  non-overlapping  segments  and  rectangular  windows  are  used.  The  variance  of  the  diagonal  elements 
is  twice  as  large  as  that  of  the  off-diagonal  elements. 

A  computer  simulation  was  performed  to  test  the  accuracy  of  the  preceding  theoretical  result.  The 
bispectrum  of  a  white,  Gaussian  noise  process  of  variance  0.25  Volts"  was  calculated.  The  sampling 
frequency  is  1/2  Hertz.  The  corresponding  power  density  spectrum  is  theoretically  flat  and  equal  to 
0.5  Watts/Hertz.  Rectangular  windows  and  non-overlapping  segments  were  used.  According  to  the  above 
equation  the  variance  of  the  off-diagonal  elements  of  the  bispectrum  is  equal  to 

var  {B,„(/i,/a)}  =  2^.  (79) 

Figure  4  shows  a  graph  of  this  equation  along  with  the  simulation  results  (shown  as  points)  as  a  function 
of  D  for  a  value  of  N  equal  to  4096.  Figure  5  graphs  the  variance  as  a  function  of  N  for  a  value  of  D 
equal  to  256.  (Only  values  of  D  and  N  that  are  equal  to  a  power  of  two  were  used  in  the  simulation.) 
There  is  good  agreement  between  the  approximate  theoretical  equation  and  the  simulation  results.  The 
variance  is  also  shown  in  the  two  figures  for  the  cases  in  which  the  data  windows  were  overlapped  by  50% 
and  75%.  To  obtain  the  theoretical  curves  for  these  two  cases  the  equation  for  variance  with  no  overlap 
was  multiplied  by  (1  —  a/2)  where  a  is  the  fractional  overlap.  (This  is  just  a  rule-of-thumb  that  seems  to 
work  well.)  As  expected  the  variance  decreases  as  the  amount  of  overlap  is  increased. 
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Figure  5:  Variance  of  the  bispectrum  of  white  Gaussian  noise  as  a  function  of  observation  length  (rectan¬ 
gular  window).  The  top,  middle  and  bottom  lines  correspond  to  0%,  50%  and  75%  overlap. 


Huber  also  gives  an  approximate  formula  for  the  variance  of  the  bicoherence  [8] 

var{6xxx(/i ,  /2)}  =  (80) 

White  Gaussian  noise  was  used  to  test  this  equation  for  the  same  set  of  parameters  used  in  the  bispectrum 
simulation.  The  results  are  shown  in  Figures  6  and  7.  When  overlapping  windows  were  used  the  same 
rule-of-thumb  used  to  adjust  the  expression  for  the  variance  of  the  bispectrum  was  applied.  Again  there 
is  very  good  agreement  between  the  formula  and  the  simulation  results.  Simulation  also  verified  that  the 
variance  of  the  bicoherence  estimate  was  relatively  independent  of  the  level  of  noise  variance. 

The  data  can  be  tapered  at  the  ends  of  the  window  to  further  reduce  the  variance.  Figs.  8  through  11 
show  simulation  results  when  a  Kaiser-Bessel  window  with  parameter  /?  =  2.5t  is  used.  These  figures 
should  be  compared  to  Figs.  4  through  7. 

Since  the  bicoherence  is  asymptotically  normal,  |6xxx(/i ,  /2)|2  will  be  asymptotically \2  distributed  with 
two  degrees  of  freedom  for  a  process  with  a  vanishing  true  bispectrum  (for  example,  a  Gaussian  process). 
Since  the  expected  value  of  \bxxx(f\ ,  /a))2  is  equal  to  var{&xxx(/i ,  /2)}  in  this  case  the  distribution  is  known 
and  we  can  determine  levels  at  which  the  magnitude  squared  bicoherence  is  statistically  significant.  For 
example,  the  probability  that  \bxxx(f\,  /2)|2  exceeds  three  times  its  mean  value  is  less  than  0.05,  i.e. 

Prob{j6xxx(/i,  /2)|2  >  3var{6xxx(/i,  /2)}}  <  0.05  (81) 

or 

Prob{|6xxx(/i,/2)|  >  >/3var{6xxx(/i,/2)}}  <  0.05.  (82) 

In  order  to  detect  quadratic  or  harmonic  coupling  it  is  necessary  for  the  significance  level  to  be  less  than 
sfUT  (see  Eqs.  75  and  77).  As  a  rule-of-thumb  we  select  the  bispectrum  estimation  algorithm  parameters 
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Figure  6:  Variance  of  the  bicoherence  of  white  Gaussian  noise  as  a  function  of  segment  length  (rectangular 
window).  The  top,  middle  and  bottom  lines  correspond  to  0%,  50%  and  75%  overlap. 
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Figure  7:  Variance  of  the  bicoherence  of  white  Gaussian  noise  as  a  function  of  observation  length  (rectan¬ 
gular  window).  The  top,  middle  and  bottom  lines  correspond  to  0%,  50%  and  75%  overlap. 
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Figure  8:  Variance  of  the  bispectrum  of  white  Gaussian  noise  as  a  function  of  segment  length  (Kaiser- 
Bessel  window).  The  top  set  of  points  (diamonds)  are  simulation  results  when  no  overlapping  was  used. 
The  middle  set  of  points  (pluses)  correspond  to  50%  overlap  and  the  bottom  points  (squares)  correspond 
to  75%  overlap. 
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Figure  9:  Variance  of  the  bispectrum  of  white  Gaussian  noise  as  a  function  of  observation  length  (Kaiser- 
Bessel  window).  The  top,  middle  and  bottom  lines  correspond  to  0%,  50%  and  75%  overlap. 
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Figure  10:  Variance  of  the  bicoherence  of  white  Gaussian  noise  as  a  function  of  segment  length  (Kaiser- 
Bessel  window).  The  top,  middle  and  bottom  lines  correspond  to  0%,  50%  and  75%  overlap. 
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Figure  11:  Variance  of  the  bicoherence  of  white  Gaussian  noise  as  a  function  of  observation  length  (Kaiser- 
Bessel  window).  The  top,  middle  and  bottom  lines  correspond  to  0%,  50%  and  75%  overlap. 


( N ,  D,  etc.)  such  that  the  variance  in  the  white  Gaussian  noise  case  satisfies 

v/3var{6  XXX  (/i,/a)}  <  v'DT/ 5 


var{6xxx(/i,/2)}  <  DT/75.  (84) 

From  Eq.  80  this  implies,  that  for  non-overlapping,  rectangular  data  windows  D  and  N  should  satisfy 

#  *  W  <85> 


For  practical  reasons  it  may  be  necessary  to  use  overlapping  and  nonrectangular  windows  in  order  to 
increase  the  ratio  of  D  to  N.  In  that  case,  simulation  is  used  to  determine  the  estimation  algorithm 
parameters  such  that  Eq.  84  is  satisfied. 
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4  Bispectral  Analysis  Results 

We  now  present  results  from  applying  bispectral  analysis  to  real  ocean  acoustic  data.  In  each  case, 
graphs  of  the  power  density  spectrum,  the  magnitude  of  the  bispectrum,  the  bispectrum  normalization 
function,  and  the  magnitude  of  the  bicoherence  are  presented.  The  bispectrum  normalization  function  is 
the  denominator  of  Eq.  11.  Graphs  of  this  function  are  useful  in  interpreting  bicoherence  results.  The 
power  density  spectrum  is  plotted  in  decibels.  The  bispectrum  and  the  bispectrum  normalization  are  both 
plotted  on  the  same  scale.  We  plot  20/3  times  the  base  10  logarithm  of  both  functions  and  the  peak  of 
the  bispectrum  normalization  function  is  defined  as  0  dB.  (The  20/3  value  is  just  a  convenient  factor  to 
use  for  comparing  bispectrum  and  power  density  spectrum  levels.)  The  magnitude  of  the  bicoherence  is 
normalized  by  dividing  by  \/DT  the  result  is  then  plotted  on  a  linear  scale  with  range  between  0  and  1. 

4.1  Detection  of  a  Non-Gaussian,  Nonlinear  Process 

The  bispectrum  of  a  data  set  from  the  August  1988  Swallow  float  experiment  [20]  was  calculated.  These 
data  were  recorded  by  float  5  which  was  at  a  depth  of  about  700  m.  The  bottom  depth  was  approximately 
3800  m.  The  Swallow  floats  contain  three  geophones  and  one  hydrophone.  Only  the  bispectrum  of  the 
hydrophone  is  presented  here.  The  Swallow  floats  use  a  sampling  rate  of  50  Hz.  A  single  record  (2250 
samples  or  45  seconds  of  data)  was  selected  for  bispectral  analysis.  To  satisfy  Eq.  84  in  the  white  Gaussian 
noise  case,  Eq.  85  suggests  a  segment  length  of  only  30  points  if  non-overlapping,  rectangular  windows  are 
used.  Simulation  indicates  that  the  variance  of  the  bicoherence  magnitudes  of  white  Gaussian  noise  with 
N  =  2240,  D  =  128,  50%  overlap  and  a  Kaiser-Bessel  window  is  equal  to  0.035.  This  satisfies  Eq.  84  and 
substituting  into  Eq.  82  gives 

Prob{|ixxx(./i,/2)|>  0.32}  <  0.05,  (86) 

or,  for  the  normalized  bicoherence  value, 

Prob{\bxxx(flth)\/VDT  >  0.2}  <  0.05  (87) 

The  results  are  shown  in  Figs.  12  through  19.  Surface  and  gray-scale  plots  are  shown  for  the  estimated 
bispectrum,  bispectrum  normalization  and  normalized  bicoherence.  Two  surface  plots  of  the  normalized 
bicoherence  are  shown.  The  second  surface  plot  has  a  raised  “floor”  at  0.2.  Points  in  the  normalized 
bicoherence  above  the  floor  are  statistically  significant  at  the  0.05  level.  In  the  normalized  bicoherence 
gray-scale  plots,  values  at  or  above  the  0.2  level  are  shown  as  white. 

There  are  approximately  1024  points  in  the  principle  domain  of  the  bicoherence.  If  the  data  were 
white  Gaussian2  we  would  expect  approximately  50  points  in  our  bicoherence  estimate  to  have  values 
above  0.2  (which  corresponds  to  the  0.05  statistical  significance  level).  There  are  actually  124  points  that 
have  values  greater  than  0.2.  We  only  expect  five  points  to  exceed  0.27  (which  corresponds  to  the  0.005 
significance  level).  There  are  58  points  in  our  estimated  bicoherence  that  exceed  this  level.  Examination 
of  the  bicoherence  estimates  shown  in  Figs.  17  through  19  indicates  that  the  low-frequency  portion  of  the 
signal  (below  5  Hz)  is  non-Gaussian  and  nonlinear.  The  Hinich  tests  [6]  could  be  applied  to  further  aid  in 
deciding  whether  or  not  the  time  series  is  in  fact  non-Gaussian  and  nonlinear. 

For  comparison  the  power  density  spectrum,  bispectrum,  and  bicoherence  of  the  next  data  record  was 
also  calculated.  The  results  are  shown  in  Figs.  20  through  27.  The  bicoherence  estimate  for  this  data 
record  looks  much  like  the  previous  one.  There  are  more  points  that  exceed  the  0.2  value  (152),  but  a 
fewer  that  exceed  the  0.27  value  (54). 

4.2  Detection  of  a  Harmonic  Process 

A  power  density  spectrum  for  a  particular  set  of  data  taken  from  the  1991  MDA  experiment  is  shown  in 
Fig.  28.  (This  particular  data  set  was  obtained  from  element  17  of  leg  A  of  the  array  on  day  196  starting  at 

2  If  the  data  were  correlated  Gaussian  we  would  expect  even  fewer  values  to  exceed  the  0.2  level. 
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Figure  12:  Power  density  spectrum  of  the  acoustic  pressure  as  measured  by  a  Swallow  float.  First  data 
record. 
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Figure  20:  Power  density  spectrum  of  the  acoustic  pressure  as  measured  by  a  Swallow  float.  Second  data 
record. 
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Figure  28:  Power  density  spectrum  of  the  acoustic  pressure  as  measured  by  element  17  of  leg  A  of  the 
MDA  array. 


time  00:52:39.483.)  This  data  set  was  selected  because  of  the  set  of  lines  that  appear  in  the  power  density 
spectrum  at  multiples  of  approximately  9.1  Hz.  Their  appearance  was  correlated  with  the  sighting  of  a 
large  ship  transiting  the  area.  The  power  density  spectrum  was  estimated  using  Welch’s  method  and  a  512 
point  FFT.  The  data  was  tapered  (using  a  Kaiser-Bessel  window  with  parameter  f3  —  2.5t)  and  averaged 
using  50%  overlap.  The  sampling  frequency  is  150  Hz.  Using  these  same  parameters,  simulation  indicates 
that  the  variance  of  the  bicoherence  of  white  Gaussian  noise  is  equal  to  0.044  for  an  observation  time  of 
64  seconds  (9600  samples).  This  satisfies  Eq.  84  and  substituting  into  Eq.  82  gives 

Prob {\bxtx(fi,f2)\  >  0.363}  <  0.05  (88) 

or 

Prob{|ixxx(/1,/2)|/v/£T>  0.2}  <  0.05.  (89) 

The  bispectrum  and  the  bispectrum  normalization  function  (the  denominator  of  the  bicoherence)  for 
this  data  set  are  shown  in  Figs.  29  through  32.  Examination  of  the  power  density  spectrum  indicates 
that  a  harmonic  process  with  a  fundamental  frequency  of  9.1  Hz  may  be  present.  There  is  a  plainly  visible 
grid  in  Fig.  30.  The  regularity  of  the  grid  spacing  provides  an  indication  that  the  process  contains  many 
harmonically  related  components. 

The  normalized  bicoherence  for  this  case  is  shown  in  Figs.  33  through  35.  In  Fig.  34  there  are  several 
peaks  above  the  0.2  bicoherence  level.  This  is  to  be  expected  because  of  the  large  number  of  points  in  the 
principle  domain  of  the  bispectrum  plane.  For  an  FFT  length  of  512  points  there  are  16512  (N(N +4)/16) 
points  in  the  principle  domain  of  the  bispectrum  plane.  For  a  white  Gaussian  noise  process  we  expect 
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Figure  29:  Bispectrum  of  the  acoustic  pressure  as  measured  by  element  17  of  leg  A  of  the  MDA  array. 
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Figure  30:  Bispectrum  of  the  acoustic  pressure  as  measured  by  element  17  of  leg  A  of  the  MDA  array. 
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Figure  34:  Bicoherence  of  the  acoustic  pressure  as  measured  by  element  17  of  leg  A  of  the  MDA  array. 
Floor  raised  to  0.2. 
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Figure  35:  Bicoherence  of  the  acoustic  pressure  as  measured  by  element  17  of  leg  A  of  the  MDA  array 


Peak 

Value 

h 

h 

1 

0.4671600 

9.0820312 

9.0820312 

2 

0.4640700 

36.6210938 

9.3750000 

3 

0.4503400 

18.4570312 

9.0820312 

4 

0.4412500 

27.5390625 

27.2460938 

5 

0.4384600 

34.8632812 

9.0820312 

6 

0.4345900 

39.8437500 

9.0820312 

7 

0.4082000 

29.8828125 

9.0820312 

8 

0.4057400 

37.5000000 

9.3750000 

9 

0.3947200 

54.4921875 

9.0820312 

10 

0.3864100 

18.1640625 

14.6484375 

11 

0.3838200 

27.2460938 

18.1640625 

12 

0.3713800 

23.7304688 

18.1640625 

13 

0.3636300 

40.1367188 

28.1250000 

14 

0.3618500 

38.3789062 

9.0820312 

15 

0.3472300 

31.0546875 

29.5898438 

16 

0.3445000 

40.4296875 

0.2929688 

17 

0.3422700 

44.2382812 

9.0820312 

18 

0.3416200 

45.9960938 

18.1640625 

19 

0.3414100 

33.9843750 

13.7695312 

20 

0.3396300 

54.7851562 

11.1328125 

Table  1:  The  20  largest  peaks  in  the  bicoherence  of  the  acoustic  pressure  as  measured  by  element  17  of 
leg  A  of  the  MDA  array. 


826  points  to  exceed  the  0.2  bicoherence  value.  There  are  actually  1273  points  in  Fig.  34  that  exceed  the 
0.2  value.  This  suggests  that  the  data  may  have  some  non-Gaussian  characteristics.  The  20  largest  peaks 
in  the  bicoherence  and  their  locations  are  listed  in  Table  1.  The  four  largest  peaks  all  occur  within  a 
single  FFT  bin  width  of  multiples  of  9.082  Hz.  The  largest  peak  has  coordinates  (9.082,  9.082)  in  the  /i, 
fi  frequency  plane.  This  is  an  indication  of  coupling  between  the  9.082  and  18.164  Hz  components.  The 
second  largest  peak  in  the  bicoherence  occurs  at  (36.621,  9.375).  The  9.375  value  is  one  FFT  bin  away  from 
the  9.082  value  and  the  36.621  value  is  one  FFT  bin  away  from  four  times  9.082.  This  indicates  coupling 
between  the  fundamental,  the  third,  and  fourth  harmonic  components.  The  third  highest  peak  occurs  at 
(18.457,  9.082)  and  indicates  coupling  between  the  fundamental,  the  first,  and  the  second  harmonics.  The 
fourth  peak  is  located  at  (27.539,  27.246)  and  implies  coupling  between  the  second  and  fifth  harmonics. 
The  next  two  peaks  occur  in  the  area  about  (36.621,  9.082)  which  is  in  the  region  of  the  third  harmonic. 
In  Fig.  34  we  can  see  considerable  spread  about  the  third  harmonic  along  the  /a  =  9.082  line.  Referring 
to  the  plot  of  the  power  density  spectrum  we  can  see  that  the  third  harmonic  is  much  weaker  than  the 
others.  Going  further  down  the  table  one  can  find  other  peaks  that  are  located  near  multiples  of  9.082  Hz. 
The  ninth  peak,  for  example,  occurs  at  (54.492,  9.082).  The  value  54.492  is  exactly  six  times  greater  than 
9.082  and  indicates  coupling  between  the  fundamental  and  the  fifth  and  sixth  harmonics.  Most  of  the 
time,  if  a  peak  does  not  have  coordinates  which  are  exact  multiples  of  9.082,  it  is  located  in  a  bin  with  a 
slightly  higher  frequency.  This  indicates  that  the  fundamental  frequency  is  slightly  higher  than  9.082  Hz. 

As  another  example,  data  collected  at  the  same  time  but  from  a  different  element  of  the  array  (element 
33)  were  analyzed  using  the  same  parameters  as  in  the  previous  case.  The  results  are  shown  in  Figs.  36 
through  43. 

There  are  1495  points  that  exceed  the  0.2  bicoherence  value  for  this  element.  Again,  if  the  process 
were  white  and  Gaussian  we  only  expect  826  points  to  have  values  above  0.2.  Looking  at  the  table  of 
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Figure  36:  Power  density  spectrum  of  the  acoustic  pressure  as  measured  by  element  33  of  leg  A  of  the 
MDA  array. 
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Figure  42:  Bicoherence  of  the  acoustic  pressure  as  measured  by  element  33  of  leg  A  of  the  MDA  array 
Floor  raised  to  0.2. 
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Figure  43:  Bicoherence  of  the  acoustic  pressure  as  measured  by  element  33  of  leg  A  of  the  MDA  array. 


Ltude 


Peak 

Value 

fi 

h 

1 

0.5725200 

9.3750000 

9.0820312 

2 

0.4942200 

18.4570312 

9.0820312 

3 

0.4583200 

19.3359375 

12.8906250 

4 

0.4269600 

31.6406250 

6.1523438 

5 

0.4207500 

54.7851562 

9.0820312 

6 

0.4202700 

0.2929688 

0.2929688 

7 

0.3956600 

27.2460938 

9.3750000 

8 

0.3853800 

45.4101562 

9.0820312 

9 

0.3849600 

19.0429688 

13.4765625 

10 

0.3826200 

20.8007812 

11.4257812 

11 

0.3709000 

48.6328125 

0.2929688 

12 

0.3685600 

13.1835938 

12.8906250 

13 

0.3659900 

34.5703125 

17.8710938 

14 

0.3571500 

36.6210938 

36.3281250 

15 

0.3545500 

53.3203125 

15.8203125 

16 

0.3530600 

18.1640625 

9.6679688 

17 

0.3527000 

13.1835938 

10.5468750 

18 

0.3519800 

57.1289062 

9.0820312 

19 

0.3518800 

37.5000000 

31.6406250 

20 

0.3511600 

16.1132812 

15.5273438 

Table  2:  The  20  largest  peaks  in  the  bicoherence  of  the  acoustic  pressure  as  measured  by  element  33  of 
leg  A  of  the  MDA  array. 
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peak  values  and  their  locations  we  can  again  see  that  several  of  the  peaks  have  coordinates  that  are  near 
multiples  of  9.082.  This  provides  further  evidence  that  component  lines  in  the  power  density  spectrum  are 
harmonically  related. 
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5  Conclusions 


Bispectral  analysis  appears  to  be  a  promising  technique  for  use  in  detecting  non-Gaussianity  and  non¬ 
linearity.  Although  the  bispectrum  of  a  harmonic  process  is  zero  the  particular  bispectrum  estimation 
algorithm  discussed  here  has  certain  properties  that  make  it  useful  in  detecting  harmonics. 

The  bicoherence  plays  a  key  role  in  the  detection  of  non-Gaussian  and  non-linear  characteristics  and 
also  in  detecting  frequency  coupling.  Knowledge  of  the  variance  of  the  estimated  bicoherence  under  the 
Gaussian  hypothesis  is  required  to  determine  if  peaks  in  the  bicoherence  are  statistically  significant.  For 
the  estimation  algorithm  discussed  here,  a  simple  expression  for  the  variance  exists  when  rectangular 
windows  are  used.  For  non-rectangular  windows,  a  simulation  is  required  to  determine  the  variance. 
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